Oyster exposure experiments, length and weight analyses
knitr::opts_chunk$set(echo = TRUE, warning = F, message = F, fig.path = 'figs/', dev.args = list(family = 'serif'))
library(tidyverse)
library(janitor)
library(readxl)
library(patchwork)
library(survival)
library(ggfortify)
source("R/funcs.R")
data(allind)
data(indbymods)
data(dirindbymods)
# consistent treatment terminology names
trts <- tibble(
shrtlab = c("7.7C", "7.7A0.2", "7.7A0.5", "8.0C", "8.0A0.2"),
lngslab = c( "7.7 Constant", "7.7 Fluctuating 0.2A", "7.7 Fluctuating 0.5A", "8.0 Constant", "8.0 Fluctuating 0.2A")
)
# week combinations
wks <- sort(unique(allind$week)) %>%
combn(2) %>%
t %>%
data.frame %>%
rename(
week1 = X1,
week2 = X2
)
grds <- allind %>%
select(trt, species, var) %>%
unique() %>%
crossing(wks) %>%
mutate(
res = purrr::pmap(list(trt, species, var, week1, week2), function(trt, species, var, week1, week2){
specsel <- species
varsel <- var
trtsel <- trt
tomod <- allind %>%
filter(species %in% specsel) %>%
filter(var %in% varsel) %>%
filter(trt %in% trtsel)
val1 <- tomod %>%
filter(week == week1) %>%
pull(val)
val2 <- tomod %>%
filter(week == week2) %>%
pull(val)
est <- t.test(val1, val2)
perc <- as.numeric(100 * diff(est$estimate) / est$estimate[2])
pout <- p_ast(est$p.value)
out <- data.frame(perc, pout)
return(out)
})
) %>%
unnest('res')Size
wgtvars <- c('Final size (cm2)')
toplo1 <- allind %>%
filter(var %in% wgtvars) %>%
group_by(week, trt, species, var) %>%
summarise(
aveval = mean(val, na.rm = T),
lowval = t.test(val)$conf.int[1],
uppval = t.test(val)$conf.int[2]
)
p1 <- ggplot(toplo1, aes(x = week, y = aveval)) +
geom_point(size = 2) +
geom_errorbar(aes(ymax = uppval, ymin = lowval), width = 1) +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 2, 4, 6)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Average values and 95% confidence intervals across time periods',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Average (+/- 95% CI)',
x = 'Week'
)
toplo2 <- grds %>%
filter(var %in% wgtvars)
p2 <- ggplot(toplo2, aes(x = week1, y = week2, fill = perc)) +
geom_tile(color = 'black') +
geom_text(aes(label = pout)) +
scale_x_continuous('Week', expand = c(0, 0), breaks = c(0, 2, 4)) +
scale_y_continuous('Week', expand = c(0, 0), breaks = c(2, 4, 6)) +
# scale_fill_distiller('% diff.', palette = 'Purples', direction = 1) +
scale_fill_gradient2(low = 'lightgreen', mid = 'white', high = 'lightblue') +
facet_grid(trt ~ species + var) +
theme(
panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank()
) +
labs(
title = 'Differences in means between experiment time points',
subtitle = 'Significant differences are indicated by text'
)
toplo3 <- allind %>%
filter(var %in% wgtvars)
p3 <- ggplot(toplo3, aes(x = val)) +
stat_ecdf(aes(colour = trt), size = 1) +
facet_grid(week ~ species + var, scales = 'free') +
scale_colour_viridis_d() +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank()
) +
labs(
title = 'Cumulative distribution plots',
subtitle = 'Panels separated by species, week, and response measure',
y = 'Cumulative frequency',
x = 'Weight (g)'
)
p1 + p2 + p3 + plot_layout(ncol = 2)
Size across individuals
Observed
# all individuals
datprp <- read_csv(here::here('data/raw/finaldata-with-newmissingdata.csv')) %>%
clean_names() %>%
select(week, trt = treatment, jar, id = individual_id, species, final_size = surface_area_cm2_final, initial_size = surface_area_cm2_initial)
toplo1 <- datprp %>%
gather('var', 'val', initial_size, final_size) %>%
separate(var, c('period', 'var')) %>%
group_by(week, trt, jar, id, species, period, var) %>%
summarise(val = mean(val, na.rm = T)) %>%
ungroup %>%
mutate(
week = case_when(
period == 'initial' ~ 0,
T ~ week
)
) %>%
unite('jar_id', jar, id)
p1 <- ggplot(toplo1, aes(x = week, y = val, group = jar_id)) +
geom_line() +
geom_point(size = 2) +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 2, 4, 6)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking size by individuals from time zero',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Size (cm2)',
x = 'Week'
)
toplo2 <- datprp %>%
mutate(
delt_size = final_size - initial_size,
init_size = 0,
) %>%
filter(week != 0) %>%
select(week, trt, jar, id, species, delt_size, init_size) %>%
gather('var', 'val', delt_size, init_size) %>%
separate(var, c('period', 'var')) %>%
group_by(week, trt, jar, id, species, period, var) %>%
summarise(val = mean(val, na.rm = T)) %>%
ungroup %>%
filter(!week == 0) %>%
mutate(
week = case_when(
period == 'init' ~ 0,
T ~ week
)
) %>%
unite('jar_id', jar, id)
p2 <- ggplot(toplo2, aes(x = week, y = val, group = jar_id)) +
geom_line() +
geom_point(size = 2) +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 2, 4, 6)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking size by individuals from time zero',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Size change from time zero (g)',
x = 'Week'
) +
geom_hline(yintercept = 0, linetype = 'dotted', colour = 'black')
toplo3 <- datprp %>%
mutate(
scld_size = (final_size - initial_size) / week + initial_size,
init_size = initial_size,
) %>%
filter(week != 0) %>%
select(week, trt, jar, id, species, scld_size, init_size) %>%
gather('var', 'val', scld_size, init_size) %>%
separate(var, c('period', 'var')) %>%
group_by(week, trt, jar, id, species, period, var) %>%
summarise(val = mean(val, na.rm = T)) %>%
ungroup %>%
filter(!week == 0) %>%
mutate(
week = case_when(
period == 'init' ~ 0L,
T ~ 1L
)
) %>%
unite('jar_id', jar, id)
p3 <- ggplot(toplo3, aes(x = week, y = val, group = jar_id)) +
geom_line(size = 0.5) +
geom_point(size = 2) +
# geom_point(aes(y = aveval), colour = 'black', size = 3) +
# geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
# geom_line(aes(y = aveval), colour = 'black') +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 1)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking scaled size by individuals from time zero, all weeks',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Scaled size (cm2/week)',
x = 'Week'
)
toplo4 <- toplo3 %>%
group_by(week, trt, species, period, var) %>%
summarize(
aveval = mean(val, na.rm = T),
uprval = t.test(val)$conf.int[2],
lwrval = t.test(val)$conf.int[1]
)
p4 <- ggplot(toplo4, aes(x = week)) +
geom_point(aes(y = aveval), colour = 'black', size = 3) +
geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
geom_line(aes(y = aveval), colour = 'black') +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 1)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking scaled size by individuals from time zero, all weeks',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Average scaled size (cm2/week) +/- 95% CI',
x = 'Week'
)
toplo5 <- datprp %>%
mutate(
scld_size = (final_size - initial_size) / week + initial_size,
init_size = initial_size,
) %>%
filter(week == 2) %>%
select(week, trt, jar, id, species, scld_size, init_size) %>%
gather('var', 'val', scld_size, init_size) %>%
separate(var, c('period', 'var')) %>%
group_by(week, trt, jar, id, species, period, var) %>%
summarise(val = mean(val, na.rm = T)) %>%
ungroup %>%
mutate(
week = case_when(
period == 'init' ~ 0L,
T ~ 1L
)
) %>%
unite('jar_id', jar, id)
p5 <- ggplot(toplo5, aes(x = week, y = val, group = jar_id)) +
geom_line(size = 0.5) +
geom_point(size = 2) +
# geom_point(aes(y = aveval), colour = 'black', size = 3) +
# geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
# geom_line(aes(y = aveval), colour = 'black') +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 1)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking scaled size by individuals from time zero, week two only',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Scaled size (cm2/week)',
x = 'Week'
)
toplo6 <- toplo5 %>%
group_by(week, trt, species, period, var) %>%
summarize(
aveval = mean(val, na.rm = T),
uprval = t.test(val)$conf.int[2],
lwrval = t.test(val)$conf.int[1]
)
p6 <- ggplot(toplo6, aes(x = week)) +
geom_point(aes(y = aveval), colour = 'black', size = 3) +
geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
geom_line(aes(y = aveval), colour = 'black') +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 1)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking scaled size by individuals from time zero, week two only',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Average Scaled size (cm2/week) +/- 95% CI',
x = 'Week'
)
toplo7 <- datprp %>%
mutate(
scld_size = (final_size - initial_size) / week + initial_size,
init_size = initial_size,
) %>%
filter(week == 4) %>%
select(week, trt, jar, id, species, scld_size, init_size) %>%
gather('var', 'val', scld_size, init_size) %>%
separate(var, c('period', 'var')) %>%
group_by(week, trt, jar, id, species, period, var) %>%
summarise(val = mean(val, na.rm = T)) %>%
ungroup %>%
mutate(
week = case_when(
period == 'init' ~ 0L,
T ~ 1L
)
) %>%
unite('jar_id', jar, id)
p7 <- ggplot(toplo7, aes(x = week, y = val, group = jar_id)) +
geom_line(size = 0.5) +
geom_point(size = 2) +
# geom_point(aes(y = aveval), colour = 'black', size = 3) +
# geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
# geom_line(aes(y = aveval), colour = 'black') +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 1)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking scaled size by individuals from time zero, week four only',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Scaled size (cm2/week)',
x = 'Week'
)
toplo8 <- toplo7 %>%
group_by(week, trt, species, period, var) %>%
summarize(
aveval = mean(val, na.rm = T),
uprval = t.test(val)$conf.int[2],
lwrval = t.test(val)$conf.int[1]
)
p8 <- ggplot(toplo8, aes(x = week)) +
geom_point(aes(y = aveval), colour = 'black', size = 3) +
geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
geom_line(aes(y = aveval), colour = 'black') +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 1)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking scaled size by individuals from time zero, week four only',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Average Scaled size (cm2/week) +/- 95% CI',
x = 'Week'
)
toplo9 <- datprp %>%
mutate(
scld_size = (final_size - initial_size) / week + initial_size,
init_size = initial_size,
) %>%
filter(week == 6) %>%
select(week, trt, jar, id, species, scld_size, init_size) %>%
gather('var', 'val', scld_size, init_size) %>%
separate(var, c('period', 'var')) %>%
group_by(week, trt, jar, id, species, period, var) %>%
summarise(val = mean(val, na.rm = T)) %>%
ungroup %>%
mutate(
week = case_when(
period == 'init' ~ 0L,
T ~ 1L
)
) %>%
unite('jar_id', jar, id)
p9 <- ggplot(toplo9, aes(x = week, y = val, group = jar_id)) +
geom_line(size = 0.5) +
geom_point(size = 2) +
# geom_point(aes(y = aveval), colour = 'black', size = 3) +
# geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
# geom_line(aes(y = aveval), colour = 'black') +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 1)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking scaled size by individuals from time zero, week six only',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Scaled size (cm2/week)',
x = 'Week'
)
toplo10 <- toplo9 %>%
group_by(week, trt, species, period, var) %>%
summarize(
aveval = mean(val, na.rm = T),
uprval = t.test(val)$conf.int[2],
lwrval = t.test(val)$conf.int[1]
)
p10 <- ggplot(toplo4, aes(x = week)) +
geom_point(aes(y = aveval), colour = 'black', size = 3) +
geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
geom_line(aes(y = aveval), colour = 'black') +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 1)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking scaled size by individuals from time zero, week six only',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Average Scaled size (cm2/week) +/- 95% CI',
x = 'Week'
)
toplo11 <- datprp %>%
mutate(
delt_size = final_size - initial_size,
init_size = 0,
) %>%
filter(week != 0) %>%
select(week, trt, jar, id, species, delt_size, init_size) %>%
gather('var', 'val', delt_size, init_size) %>%
separate(var, c('period', 'var')) %>%
group_by(week, trt, jar, id, species, period, var) %>%
summarise(val = mean(val, na.rm = T)) %>%
ungroup %>%
filter(!week == 0) %>%
mutate(
week = case_when(
period == 'init' ~ 0,
T ~ week
)
) %>%
unite('jar_id', jar, id) %>%
filter(week != 0) %>%
mutate(
val = case_when(
sign(val) %in% c(0, 1) ~ 1,
sign(val) == -1 ~ -1
)
) %>%
group_by(week, trt, species, period, var, val) %>%
summarise(cnt = sum(val)) %>%
na.omit() %>%
ungroup %>%
mutate(val = factor(val, levels = c(1, -1), labels = c('positive', 'negative')))
p11 <- ggplot(toplo11, aes(x = week, y = cnt, fill = factor(val))) +
geom_bar(stat = 'identity') +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(2, 4, 6)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank()
) +
labs(
title = 'Tracking size by individuals from time zero',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Count of individuals with positive or negative size change',
x = 'Week'
)
p1 + p2 + p11 + plot_spacer() + p3 + p4 + p5 + p6 + p7 + p8 + p9 + p10 + plot_layout(ncol = 2)
Modelled, all changes
Olympia oyster
indbymods %>%
filter(species == 'Olympia') %>%
filter(var == 'Final size (cm2)') %>%
pull(plomod)
Pacific oyster
indbymods %>%
filter(species == 'Pacific') %>%
filter(var == 'Final size (cm2)') %>%
pull(plomod)
Modelled, positive size
Olympia oyster
dirindbymods %>%
filter(species == 'Olympia') %>%
filter(var == 'Final size (cm2)') %>%
filter(direction == 1) %>%
pull(plomod)
Pacific oyster
dirindbymods %>%
filter(species == 'Pacific') %>%
filter(var == 'Final size (cm2)') %>%
filter(direction == 1) %>%
pull(plomod)
Modelled, negative size
Olympia oyster
dirindbymods %>%
filter(species == 'Olympia') %>%
filter(var == 'Final size (cm2)') %>%
filter(direction == -1) %>%
pull(plomod)
Pacific oyster
dirindbymods %>%
filter(species == 'Pacific') %>%
filter(var == 'Final size (cm2)') %>%
filter(direction == -1) %>%
pull(plomod)
Weight
wgtvars <- c('Shell weight (g)', 'Tissue weight (g)')
toplo1 <- allind %>%
filter(var %in% wgtvars) %>%
group_by(week, trt, species, var) %>%
summarise(
aveval = mean(val, na.rm = T),
lowval = t.test(val)$conf.int[1],
uppval = t.test(val)$conf.int[2]
)
p1 <- ggplot(toplo1, aes(x = week, y = aveval)) +
geom_point(size = 2) +
geom_errorbar(aes(ymax = uppval, ymin = lowval), width = 1) +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 2, 4, 6)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Average values and 95% confidence intervals across time periods',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Average (+/- 95% CI)',
x = 'Week'
)
toplo2 <- grds %>%
filter(var %in% wgtvars)
p2 <- ggplot(toplo2, aes(x = week1, y = week2, fill = perc)) +
geom_tile(color = 'black') +
geom_text(aes(label = pout)) +
scale_x_continuous('Week', expand = c(0, 0), breaks = c(0, 2, 4)) +
scale_y_continuous('Week', expand = c(0, 0), breaks = c(2, 4, 6)) +
# scale_fill_distiller('% diff.', palette = 'Purples', direction = 1) +
scale_fill_gradient2(low = 'lightgreen', mid = 'white', high = 'lightblue') +
facet_grid(trt ~ species + var) +
theme(
panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank()
) +
labs(
title = 'Differences in means between experiment time points',
subtitle = 'Significant differences are indicated by text'
)
toplo3 <- allind %>%
filter(var %in% wgtvars)
p3 <- ggplot(toplo3, aes(x = val)) +
stat_ecdf(aes(colour = trt), size = 1) +
facet_grid(week ~ species + var, scales = 'free') +
scale_colour_viridis_d() +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank()
) +
labs(
title = 'Cumulative distribution plots',
subtitle = 'Panels separated by species, week, and response measure',
y = 'Cumulative frequency',
x = 'Weight (g)'
)
p1 + p2 + p3 + plot_layout(ncol = 2)
Weight across individuals
Observed
# calc initial final from sum of tissue and shell weight
datprp <- read.csv(here::here('data/raw/Weight_with dead but not multiple_Kelp_4_3.csv'), stringsAsFactors = F) %>%
clean_names() %>%
# filter(week != 0) %>%
filter(species != '') %>%
mutate(
species = gsub('\\s*$', '', species),
whole_organism_weight = case_when(
is.na(whole_organism_weight) & !is.na(shell_weight) & !is.na(tissue_weight) ~ shell_weight + tissue_weight,
T ~ whole_organism_weight
)
) %>%
filter(!is.na(initial_wetweight) & !is.na(whole_organism_weight)) %>%
select(week, trt = treatment, jar = i_jar, id = individual_id, species, initial_weight = initial_wetweight,
final_weight = whole_organism_weight)
toplo1 <- datprp %>%
gather('var', 'val', initial_weight, final_weight) %>%
separate(var, c('period', 'var')) %>%
group_by(week, trt, jar, id, species, period, var) %>%
summarise(val = mean(val, na.rm = T)) %>%
ungroup %>%
mutate(
week = case_when(
period == 'initial' ~ 0L,
T ~ week
)
) %>%
unite('jar_id', jar, id)
p1 <- ggplot(toplo1, aes(x = week, y = val, group = jar_id)) +
geom_line() +
geom_point(size = 2) +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 2, 4, 6)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking weight by individuals from time zero',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Weight (g)',
x = 'Week'
)
toplo2 <- datprp %>%
mutate(
delt_weight = final_weight - initial_weight,
init_weight = 0,
) %>%
filter(week != 0) %>%
select(week, trt, jar, id, species, delt_weight, init_weight) %>%
gather('var', 'val', delt_weight, init_weight) %>%
separate(var, c('period', 'var')) %>%
group_by(week, trt, jar, id, species, period, var) %>%
summarise(val = mean(val, na.rm = T)) %>%
ungroup %>%
filter(!week == 0) %>%
mutate(
week = case_when(
period == 'init' ~ 0L,
T ~ week
)
) %>%
unite('jar_id', jar, id)
p2 <- ggplot(toplo2, aes(x = week, y = val, group = jar_id)) +
geom_line() +
geom_point(size = 2) +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 2, 4, 6)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking weight by individuals from time zero',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Weight change from time zero (g)',
x = 'Week'
) +
geom_hline(yintercept = 0, linetype = 'dotted', colour = 'black')
toplo3 <- datprp %>%
mutate(
scld_weight = (final_weight - initial_weight) / week + initial_weight,
init_weight = initial_weight,
) %>%
filter(week != 0) %>%
select(week, trt, jar, id, species, scld_weight, init_weight) %>%
gather('var', 'val', scld_weight, init_weight) %>%
separate(var, c('period', 'var')) %>%
group_by(week, trt, jar, id, species, period, var) %>%
summarise(val = mean(val, na.rm = T)) %>%
ungroup %>%
filter(!week == 0) %>%
mutate(
week = case_when(
period == 'init' ~ 0L,
T ~ 1L
)
) %>%
unite('jar_id', jar, id)
p3 <- ggplot(toplo3, aes(x = week, y = val, group = jar_id)) +
geom_line(size = 0.5) +
geom_point(size = 2) +
# geom_point(aes(y = aveval), colour = 'black', size = 3) +
# geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
# geom_line(aes(y = aveval), colour = 'black') +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 1)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking scaled weight by individuals from time zero, all weeks',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Scaled weight (g/week)',
x = 'Week'
)
toplo4 <- toplo3 %>%
group_by(week, trt, species, period, var) %>%
summarize(
aveval = mean(val, na.rm = T),
uprval = t.test(val)$conf.int[2],
lwrval = t.test(val)$conf.int[1]
)
p4 <- ggplot(toplo4, aes(x = week)) +
geom_point(aes(y = aveval), colour = 'black', size = 3) +
geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
geom_line(aes(y = aveval), colour = 'black') +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 1)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking scaled weight by individuals from time zero, all weeks',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Average scaled weight (g/week) +/- 95% CI',
x = 'Week'
)
toplo5 <- datprp %>%
mutate(
scld_weight = (final_weight - initial_weight) / week + initial_weight,
init_weight = initial_weight,
) %>%
filter(week == 2) %>%
select(week, trt, jar, id, species, scld_weight, init_weight) %>%
gather('var', 'val', scld_weight, init_weight) %>%
separate(var, c('period', 'var')) %>%
group_by(week, trt, jar, id, species, period, var) %>%
summarise(val = mean(val, na.rm = T)) %>%
ungroup %>%
mutate(
week = case_when(
period == 'init' ~ 0L,
T ~ 1L
)
) %>%
unite('jar_id', jar, id)
p5 <- ggplot(toplo5, aes(x = week, y = val, group = jar_id)) +
geom_line(size = 0.5) +
geom_point(size = 2) +
# geom_point(aes(y = aveval), colour = 'black', size = 3) +
# geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
# geom_line(aes(y = aveval), colour = 'black') +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 1)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking scaled weight by individuals from time zero, week two only',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Scaled weight (g/week)',
x = 'Week'
)
toplo6 <- toplo5 %>%
group_by(week, trt, species, period, var) %>%
summarize(
aveval = mean(val, na.rm = T),
uprval = t.test(val)$conf.int[2],
lwrval = t.test(val)$conf.int[1]
)
p6 <- ggplot(toplo6, aes(x = week)) +
geom_point(aes(y = aveval), colour = 'black', size = 3) +
geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
geom_line(aes(y = aveval), colour = 'black') +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 1)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking scaled weight by individuals from time zero, week two only',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Average Scaled weight (g/week) +/- 95% CI',
x = 'Week'
)
toplo7 <- datprp %>%
mutate(
scld_weight = (final_weight - initial_weight) / week + initial_weight,
init_weight = initial_weight,
) %>%
filter(week == 4) %>%
select(week, trt, jar, id, species, scld_weight, init_weight) %>%
gather('var', 'val', scld_weight, init_weight) %>%
separate(var, c('period', 'var')) %>%
group_by(week, trt, jar, id, species, period, var) %>%
summarise(val = mean(val, na.rm = T)) %>%
ungroup %>%
mutate(
week = case_when(
period == 'init' ~ 0L,
T ~ 1L
)
) %>%
unite('jar_id', jar, id)
p7 <- ggplot(toplo7, aes(x = week, y = val, group = jar_id)) +
geom_line(size = 0.5) +
geom_point(size = 2) +
# geom_point(aes(y = aveval), colour = 'black', size = 3) +
# geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
# geom_line(aes(y = aveval), colour = 'black') +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 1)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking scaled weight by individuals from time zero, week four only',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Scaled weight (g/week)',
x = 'Week'
)
toplo8 <- toplo7 %>%
group_by(week, trt, species, period, var) %>%
summarize(
aveval = mean(val, na.rm = T),
uprval = t.test(val)$conf.int[2],
lwrval = t.test(val)$conf.int[1]
)
p8 <- ggplot(toplo8, aes(x = week)) +
geom_point(aes(y = aveval), colour = 'black', size = 3) +
geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
geom_line(aes(y = aveval), colour = 'black') +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 1)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking scaled weight by individuals from time zero, week four only',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Average Scaled weight (g/week) +/- 95% CI',
x = 'Week'
)
toplo9 <- datprp %>%
mutate(
scld_weight = (final_weight - initial_weight) / week + initial_weight,
init_weight = initial_weight,
) %>%
filter(week == 6) %>%
select(week, trt, jar, id, species, scld_weight, init_weight) %>%
gather('var', 'val', scld_weight, init_weight) %>%
separate(var, c('period', 'var')) %>%
group_by(week, trt, jar, id, species, period, var) %>%
summarise(val = mean(val, na.rm = T)) %>%
ungroup %>%
mutate(
week = case_when(
period == 'init' ~ 0L,
T ~ 1L
)
) %>%
unite('jar_id', jar, id)
p9 <- ggplot(toplo9, aes(x = week, y = val, group = jar_id)) +
geom_line(size = 0.5) +
geom_point(size = 2) +
# geom_point(aes(y = aveval), colour = 'black', size = 3) +
# geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
# geom_line(aes(y = aveval), colour = 'black') +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 1)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking scaled weight by individuals from time zero, week six only',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Scaled weight (g/week)',
x = 'Week'
)
toplo10 <- toplo9 %>%
group_by(week, trt, species, period, var) %>%
summarize(
aveval = mean(val, na.rm = T),
uprval = t.test(val)$conf.int[2],
lwrval = t.test(val)$conf.int[1]
)
p10 <- ggplot(toplo4, aes(x = week)) +
geom_point(aes(y = aveval), colour = 'black', size = 3) +
geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
geom_line(aes(y = aveval), colour = 'black') +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(0, 1)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank()
) +
labs(
title = 'Tracking scaled weight by individuals from time zero, week six only',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Average Scaled weight (g/week) +/- 95% CI',
x = 'Week'
)
toplo11 <- datprp %>%
mutate(
delt_weight = final_weight - initial_weight,
init_weight = 0,
) %>%
filter(week != 0) %>%
select(week, trt, jar, id, species, delt_weight, init_weight) %>%
gather('var', 'val', delt_weight, init_weight) %>%
separate(var, c('period', 'var')) %>%
group_by(week, trt, jar, id, species, period, var) %>%
summarise(val = mean(val, na.rm = T)) %>%
ungroup %>%
filter(!week == 0) %>%
mutate(
week = case_when(
period == 'init' ~ 0L,
T ~ week
)
) %>%
unite('jar_id', jar, id) %>%
filter(week != 0) %>%
mutate(
val = case_when(
sign(val) %in% c(0, 1) ~ 1,
sign(val) == -1 ~ -1
)
) %>%
group_by(week, trt, species, period, var, val) %>%
summarise(cnt = sum(val)) %>%
na.omit() %>%
ungroup %>%
mutate(val = factor(val, levels = c(1, -1), labels = c('positive', 'negative')))
p11 <- ggplot(toplo11, aes(x = week, y = cnt, fill = factor(val))) +
geom_bar(stat = 'identity') +
facet_grid(trt ~ species + var) +
scale_x_continuous(breaks = c(2, 4, 6)) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank()
) +
labs(
title = 'Tracking weight by individuals from time zero',
subtitle = 'Panels separated by species, treatment, and response measure',
y = 'Count of individuals with positive or negative weight change',
x = 'Week'
)
p1 + p2 + p11 + plot_spacer() + p3 + p4 + p5 + p6 + p7 + p8 + p9 + p10 + plot_layout(ncol = 2)
Modelled, all changes
Olympia oyster
indbymods %>%
filter(species == 'Olympia') %>%
filter(var == 'Whole weight (g)') %>%
pull(plomod)
Pacific oyster
indbymods %>%
filter(species == 'Pacific') %>%
filter(var == 'Whole weight (g)') %>%
pull(plomod)
Modelled, positive weight
Olympia oyster
dirindbymods %>%
filter(species == 'Olympia') %>%
filter(var == 'Whole weight (g)') %>%
filter(direction == 1) %>%
pull(plomod)
Pacific oyster
dirindbymods %>%
filter(species == 'Pacific') %>%
filter(var == 'Whole weight (g)') %>%
filter(direction == 1) %>%
pull(plomod)
Modelled, negative weight
Olympia oyster
dirindbymods %>%
filter(species == 'Olympia') %>%
filter(var == 'Whole weight (g)') %>%
filter(direction == -1) %>%
pull(plomod)
Pacific oyster
dirindbymods %>%
filter(species == 'Pacific') %>%
filter(var == 'Whole weight (g)') %>%
filter(direction == -1) %>%
pull(plomod)
Size vs weight vs. length
toplo1 <- allind %>%
group_by_at(vars(-val)) %>%
summarise(val = mean(val, na.rm = T)) %>%
ungroup %>%
spread(var, val)
p1 <- ggplot(toplo1, aes(x = `Final size (cm2)`, y = `Shell weight (g)`, colour = trt)) +
geom_point(alpha = 0.7) +
# geom_smooth(se = F, method = 'lm') +
scale_colour_viridis_d() +
facet_grid(week ~ species) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank()
) +
labs(
title = 'Tracking shell weight vs size',
subtitle = 'Panels separated by species, week, and treatment',
y = 'Shell weight (g)',
x = 'Final size (cm2)'
)
p2 <- ggplot(toplo1, aes(x = `Final size (cm2)`, y = `Tissue weight (g)`, colour = trt)) +
geom_point(alpha = 0.7) +
# geom_smooth(se = F, method = 'lm') +
scale_colour_viridis_d() +
facet_grid(week ~ species) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank()
) +
labs(
title = 'Tracking tissue weight by size',
subtitle = 'Panels separated by species, week, and treatment',
y = 'Tissue weight (g)',
x = 'Final size (cm2)'
)
p3 <- ggplot(toplo1, aes(x = `Shell weight (g)`, y = `Tissue weight (g)`, colour = trt)) +
geom_point(alpha = 0.7) +
# geom_smooth(se = F, method = 'lm') +
scale_colour_viridis_d() +
facet_grid(week ~ species) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank()
) +
labs(
title = 'Tracking tissue weight by shell weight',
subtitle = 'Panels separated by species, week, and treatment',
y = 'Tissue weight (g)',
x = 'Shell weight (g)'
)
lendat <- alllen <- read_excel(here::here('data/raw/Length_kelp_4_3_2020.xlsx')) %>%
clean_names() %>%
select(week, trt = treatment, jar, id = individual_id, species, final_length = length_cm_final, initial_length = length_cm_initial, final_width = width_cm_final, initial_width = width_cm_initial) %>%
mutate(
final_length = case_when(
week == 0 & is.na(final_length) ~ initial_length,
T ~ final_length
),
final_width = case_when(
week == 0 & is.na(final_width) ~ initial_width,
T ~ final_width
),
delt_length = final_length - initial_length,
delt_width = final_width - initial_width,
rate_length = case_when(
week == 2 ~ delt_length / 14,
week == 4 ~ delt_length / 28
),
rate_width = case_when(
week == 2 ~ delt_width / 14,
week == 4 ~ delt_width / 28
)
) %>%
gather('var', 'val', -week, -trt, -jar, -id, -species) %>%
filter(week != 6) %>%
filter(!(week %in% 0 & var %in% c('delt_length', 'delt_width', 'rate_width', 'rate_length'))) %>%
filter(!var %in% c('initial_length', 'initial_width')) %>%
group_by_at(vars(-val)) %>%
summarise(val = mean(val, na.rm = T)) %>%
ungroup %>%
spread(var, val) %>%
mutate(jar = as.character(jar))
toplo2 <- allind %>%
group_by_at(vars(-val)) %>%
summarise(val = mean(val, na.rm = T)) %>%
ungroup %>%
spread(var, val) %>%
mutate(
jar = as.character(jar),
trt = as.character(trt)
) %>%
inner_join(lendat, ., by = c('week', 'trt', 'jar', 'id', 'species')) %>%
mutate(
trt = factor(trt, levels = trts$shrtlab, labels = trts$shrtlab)
)
p4 <- ggplot(toplo2, aes(x = `Final size (cm2)`, y = final_length, colour = trt)) +
geom_point(alpha = 0.7) +
# geom_smooth(se = F, method = 'lm') +
scale_colour_viridis_d() +
facet_grid(week ~ species) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank()
) +
labs(
title = 'Tracking length by size',
subtitle = 'Panels separated by species, week, and treatment',
y = 'Final length (cm)',
x = 'Final size (cm2)'
)
p5 <- ggplot(toplo2, aes(x = `Final size (cm2)`, y = final_width, colour = trt)) +
geom_point(alpha = 0.7) +
# geom_smooth(se = F, method = 'lm') +
scale_colour_viridis_d() +
facet_grid(week ~ species) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank()
) +
labs(
title = 'Tracking width by size',
subtitle = 'Panels separated by species, week, and treatment',
y = 'Final width (cm)',
x = 'Final size (cm2)'
)
p6 <- ggplot(toplo2, aes(x = final_length, y = final_width, colour = trt)) +
geom_point(alpha = 0.7) +
# geom_smooth(se = F, method = 'lm') +
scale_colour_viridis_d() +
facet_grid(week ~ species) +
theme_bw() +
theme(
# panel.background = element_blank(),
strip.background = element_blank(),
axis.ticks = element_blank(),
panel.grid.minor = element_blank(),
legend.title = element_blank()
) +
labs(
title = 'Tracking width by length',
subtitle = 'Panels separated by species, week, and treatment',
y = 'Final width (cm)',
x = 'Final length (cm)'
)
p1 + p2 + p3 + p4 + p5 + p6 + plot_layout(ncol = 2)
Survival
srvdat <- read.csv(here::here('data/raw/Weight_with dead but not multiple_Kelp_4_3.csv'), stringsAsFactors = F) %>%
clean_names() %>%
filter(species != '') %>%
mutate(species = gsub('\\s*$', '', species)) %>%
rename(
id = individual_id,
jar = i_jar,
trt = treatment
) %>%
select(week, trt, jar, id, species, dead) %>%
mutate(
dead = case_when(
dead %in% c('x', 'X') ~ 1,
T ~ 0
),
trt = factor(trt, levels = trts$shrtlab)
) %>%
group_by(species) %>%
nest %>%
mutate(
srvmod = purrr::map(data, function(x){
fit <- survfit(Surv(week, dead) ~ trt, data = x)
return(fit)
}),
srvdif = purrr::map(data, function(x){
fit <- survdiff(Surv(week, dead) ~ trt, data = x)
return(fit)
})
)
srvplos <- srvdat %>%
mutate(
plodat = purrr::map(srvmod, fortify)
) %>%
select(-data, -srvmod, -srvdif) %>%
unnest(plodat)
cols <- c('khaki4', 'khaki3', 'khaki1', 'seagreen4', 'seagreen2')
names(cols) <- trts$shrtlab
p <- ggplot(srvplos, aes(x = time, y = surv)) +
geom_line(aes(colour = strata), size = 1.5) +
geom_ribbon(aes(ymin = lower, ymax = upper), colour = NA, fill = 'lightgrey', alpha = 0.25) +
facet_grid(species ~ strata) +
scale_colour_manual('Treatment', values = cols, guide = F) +
theme_minimal() +
theme(
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
legend.position = 'top'
) +
labs(
x = 'Week',
y = '% survival',
title = 'Survival estimates by treatment and species',
subtitle = 'Shaded region is 95% CI'
)
p
Olympia survival difference test
srvdat %>%
filter(species %in% 'Olympia') %>%
pull(srvdif) %>%
.[[1]]## Call:
## survdiff(formula = Surv(week, dead) ~ trt, data = x)
##
## n=1624, 6 observations deleted due to missingness.
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## trt=7.7C 323 30 27.0 0.3250 0.4300
## trt=7.7A0.2 324 21 27.0 1.3499 1.7860
## trt=7.7A0.5 329 29 26.8 0.1754 0.2316
## trt=8.0C 320 26 26.6 0.0125 0.0165
## trt=8.0A0.2 328 29 27.5 0.0800 0.1064
##
## Chisq= 2.1 on 4 degrees of freedom, p= 0.7Pacific survival difference test
srvdat %>%
filter(species %in% 'Pacific') %>%
pull(srvdif) %>%
.[[1]]## Call:
## survdiff(formula = Surv(week, dead) ~ trt, data = x)
##
## n=1540, 5 observations deleted due to missingness.
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## trt=7.7C 303 22 20.9 0.0595 0.0770
## trt=7.7A0.2 301 21 22.0 0.0455 0.0596
## trt=7.7A0.5 319 21 22.6 0.1136 0.1498
## trt=8.0C 325 19 22.9 0.6541 0.8647
## trt=8.0A0.2 292 27 21.6 1.3253 1.7286
##
## Chisq= 2.3 on 4 degrees of freedom, p= 0.7